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We present a new method for investigating first-order pliase transitions using Monte Carlo simu- 
lations. It relies on the multiple-histogram method and uses solely histograms of individual phases. 
In addition, we extend the method to include histograms of subphases. 

The free energy difference between phases, necessary for attributing the correct statistical weights 
to the histograms, is determined by a detour in control parameter space via auxiliary systems with 
short relaxation times. We apply this method to a recently introduced model for structure formation 
in polypeptides for which other methods fail. 
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^ I At first-order phase transitions, Monte Carlo (MC) 
^ O ' simulations encounter a particular problem of critical 
[ slowing down: since transitions between the co-existing 

■ phases are very rare, the relaxation time is extremely 
' long and usually increases exponentially with system size. 

In order to obtain equilibrated data via MC simulations 
for a system of appreciable size near the transition, pro- 
^ , hibitively long simulations would have to be performed. 

A convenient way to visualize the problem as well as 
to analyze simulation results is to employ histograms. 
Since the Hamiltonian and most observables of inter- 
est are functions of one or a few order parameters, such 
as magnetization or internal energy, a histogram of the 
0^ ' frequency of their occurence, recorded during the sim- 
""^Jj, ulation, is sufRcient to generate all information of inter- 
• est. Using the single- and the multiple-histogram method 
O ; (SHM/MHM) §, one or several histograms determined 
^ • for a system at specific sets of control parameters can 

■ be used efficiently to predict the system behaviour over 
'"O ■ a wide range of control parameter values. In such his- 

■ tograms the co-existence of phases corresponds to several 
peaks, see e.g. Fig. |^. The infrequent switching between 
individual phases makes it difficult to generate a single 

. , equilibrated histogram covering all phases, i.e. showing 

■ all peaks correctly. 
Most approaches to this problem introduce suitable 

changes to the system Hamiltonian so that transition 
states, which are encountered only rarely in the origi- 
nal system, will be sampled much more frequently in a 
simulation, thereby leading also to a higher frequency of 
switching between phases Common to all these 

methods is the idea that in a single simulation a his- 
togram can be determined which covers all interesting 
regions of state space. A drawback is that possibly 
a plethora of parameters, which render the necessary 
changes in the system Hamiltonian, have to be optimized. 

In this contribution we want to present a new approach 
to the problem. It is based on the realization that an 
equilibrated histogram of a system confined to an individ- 
ual phase, i.e. an individual peak in the full histogram, 
can be generated much more easily. Under those con- 
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ditions the relaxation times are usually small, at least 
much smaller than the switching times between phases 
in the critical regime. For a first-order phase transition, 
these single-phase histograms already cover all relevant 
system states because, even at the critical point, tran- 
sition states between the phases are encountered only 
very infrequently, and, hence, their contribution to the 
partition function is negligible. Several single-phase his- 
tograms can, therefore, be combined by the MHM to re- 
construct the complete histograms at or near the first- 
order transition. 

However, there is a technical problem involved: the 
relative contributions of various histograms, in effect the 
free energy differences between them, have to be deter- 
mined. This is a notoriously difficult task. A prerequisite 
for their correct determination, e.g. by the acceptance 
ratio method (ARM) developed by Bennett in his defini- 
tive treatment of the subject is that the histograms 
in question overlap at least partially. This is usually not 
the case, particularly not for systems with strong first- 
order transitions. To overcome this difficulty, we use a 
detour in control parameter space via auxiliary systems 
with much shorter relaxation times, which resemble the 
actual system to an extent sufficient to guarantee signif- 
icant overlaps between actual and auxiliary histograms 

We will apply this scheme to the temperature-driven 
first-order phase transition in a recently introduced sim- 
ple model of secondary and tertiary structure formation 
in polypeptides § , where the methods of Ref. fail. 

The conformation of a polypeptide of length L is rep- 
resented by a string of local conformations tTi = h, c'^, or 
c°, i = 1, . . . , L. The local conformation h corresponds to 
residues with dihedral angles characteristic of a helices. 
Any three successive monomers in helical conformation 
are spanned by a hydrogen bond with energy Ehb < 0. 
A string of / consecutive /i-residues {I > 3) forms a he- 
lix of length I. The c" and c"*" residues denote random 
coil conformations. Two helices separated solely by 
residues are taken to interact with each other, whereas 
helices with at least one residue with conformation other 
than c° between them are not. The number of contacts 
between two interacting helices, which determines the in- 
teraction energy, is taken to be equal to the length of the 
shorter helix. We set the interaction energy parameter 
to fc = 0.6 Ehb- Since the conformation space volume 
il{h) accessible to h residues is smaller than that for 
non-h residues, the conversion c"*" <-> h, for example, 
is accompanied by a change in conformational entropy. 
We define AS{(Ti) — fc^ ln[ri(cri)/r2(c+)], and assume 
equal statistical weights for the two random coil confor- 
mations, i.e., AS{c°) = AS{c+) = 0. The helix-coil 
transition in experimentally studied homo (poly) amino 
acids is adequately described by a value of AS{h)/kB = 
—4.26 -I- ln(2) « —3.57, derived from experimental data 
first discussed by Zimm and Bragg 

Using the order parameter vector of a chain conforma- 
tion {a,}, S{{a,}) = iNHBi{cT^}), Ni,{{a,}), N,i{a,})), 
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where Nhb, Nh, and Nc denote the number of hy- 
drogen bonds, /i-residues, and contacts, respectively, 
and the corresponding control parameter vector — 
{PEhb, -AS{h), f3k), with (3 = {kBT)~^, we obtain 
a free energy of the conformation {ai}, 

pF{{a,})^K-Si{a,}). (1) 

The role of F is equivalent to that of a Hamiltonian in 
other models. 

This model has a coil-dominated high-temperature 
phase and a helix-dominated low-temperature phase. For 
fc > the low-temperature phase is a single helix span- 
ning the whole system, and the change from coil to helix 
is not accompanied by a phase transition For 
fc < 0, however, the low-temperature phase is multi- 
helical, with neighboring helices stabilizing each other 
via tertiary interactions, and the transition between the 
coil and the multi-helical phases is of first order. 

The multi-helical phase consists of several subphases 
characterized by different numbers of helices. These sub- 
phases are separated from each other, and from the coil 
phase, by significant barriers in a feature which leads 
to long relaxation times at the transition point as well as 
to a slow, glass-like, relaxation within the multi-helical 
phase 1^. The existence of these subphases indicates 
that, in addition to the order parameters that enter the 
Hamiltonian, Eq. (|^), the number of helices of a confor- 
mation {(Ti}, Nhei{{o'i}), is an additional relevant order 
parameter, see also Fig. 0. 

In order to investigate the behavior of this system at 
and around its phase transition point we generated or- 
der parameter histograms, fT.i<-(5), from MC simulations 
for parameter sets corresponding to the coil and the 
multi- helical phase for various values of L. For L = 40 
additional simulations close to the phase transition point 
were performed for a check of the method. We provide 
the respective parameters and histogram properties for 
that length in TableJ. 

Using the MHM the relative probability to find the 
system at parameter set if in a state with order param- 
eter S_ , Pk{S_), can be approximated from histograms 
obtained at the parameter sets K_^, i — 1, . . . , M, by 



= ^ . .1 ■ (2) 



EtLi 9-'N, exp [-K^ ■ S + fiK^)] 



Here, Ni = ni{S_), gi — 1 + Iti renormalizes the size 
of the histogram to take correlation effects into account, 
with Ti being the correlation time of the simulation i. By 
we denote the free energy of the system at param- 
eter set K, exp[-/(^)] = Z(K) = Y.sPk(S), which is 
determined only up to an additive constant. 

We have argued above that, using Eq. (^), histograms 
of the coil and the multi-helical phase (histograms 1 and 
3 of Tab. H) suffice to reconstruct the first order phase 
transition. However, this requires the knowledge of the 
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free energy difference between these histograms. Apply- 
ing Bennett's ARM ^ to the model studied here leads 
to the equation 



{Hs-[K^-m-c))^ 



^fv = In r~' ^'x + C, (3) 



where A/y = f{K,) - f{K^), and C = HZ,N,/Z,Nj) 
has to be determined self-consistently. !F{x) — 1/[1 + 
exp(a:)] is the Fermi function, and (x) ■ denotes the aver- 
age of X with respect to histogram i. In order to apply 
Eq. (|^), the two histograms in question have to overlap 
at least partially, which is not the case for the single- 
phase histograms 1 and 3. Bennett has already pointed 
out the possibility to obtain free energy differences for 
disjoint histograms by performing additional simulations 
for intermediate parameter sets, so as to form an over- 
lapping chain of histograms. However, a simulation near 
the critical point is not feasible in general. 

In our model we can exploit the property that, by 
changing the parameter AS{h) to less negative values, 
systems are obtained which exhibit significantly shorter 
relaxation times and a more gradual transition. For such 
systems, equilibrated histograms can be generated for the 
separate phases and in the transition region (see Tab. |) 
at a much smaller expense of computation time than for 
the original system. Using such auxiliary histograms, 
a sequence of mutually overlapping histograms can be 
formed to join the single-phase histograms of the original 
system. This allows the determination of A/31 between 
the disjoint histograms 1 and 3 of Tab. | by repeated 
application of Eq. (||), as illustrated in detail in Tab. 
0. For comparison, the value for A/31 determined using 
histogram 2 of Tab. | as intermediate histogram is also 
given in Tab. ||. Both values agree closely. However, the 
expected error for the latter value is significantly larger 
due to the long relaxation time for the simulation near 
the critical point. 

It was already noted that the multi-helical phase is 
also characterized by long relaxation times, here due to 
infrequent switches between various coexisting subphases 
characterized by different helix numbers. These relax- 
ation times become prohibitively long particularly for 
systems with larger sizes (we have studied systems up 
to L = 200). With minor modifications, it is possible 
to apply the principles we used to reconstruct the his- 
togram at the phase transition also for a reconstruction 
of histograms within the multi-helical phase. Histograms 
for helical subphases with a constraint in helix number 
Nfiei can be generated with much smaller computational 
effort than for the full histogram, see Tab. |. They are 
obtained simply by introducing an infinite energy barrier 
for all steps that attempt to change the helix number 
during the simulations. 

In direct analogy to the multiple-histogram equation, 
Eq. (H), it is possible to combine several subphase his- 
tograms nf{S_) with identical helix number Nhei = a, ob- 
tained at various parameter sets K_^, i = 1, . . . , M. The 



4 



relative probability within this subphase of a system state 
S_ (with Nhei = a), at parameter set I£, is approximated 
by 

pa(g^ _ j:fU9fr'nf{S)exp[~K-S] 

Here, the sum in the denominator ranges only over sub- 
phase histograms with Nh^i — a. We note that the free 
energies 

riKj) in Eq. (|) have now to be calculated with 
respect to subphase histograms with Nhei = a only. 

These subphase histograms can now be combined, 
again using ideas of the MHM. Following the arguments 
in Rcf. [0, and using the property that the subphases 
fully partition the order parameter space into disjoint 
patches, one arrives at 

a 

for the combined histogram. Given a full histogram of 
an auxiliary system at some convenient parameter value 
K ' , one obtains for the relative contribution of subphase 
a,p^in Eq. d), 

p]^ = pheMnK)~r{K')], (6) 

where p'^, is given by p^, = Es."k'(^)/ J2s_^k'{S). 

Using the methods described we have investigated the 
first order transition of our model. Contrary to other ap- 
proaches ||T^, we choose the fluctuations of the specific 
energy e= {NHB-EHB+Nc-k)/L, i.e. 62 ((e- (e))^), 
as the observable to monitor the transition [Ol . For first 
order transitions the peak of this function assumes a fi- 
nite nonzero value in the thermodynamic limit, which is 
given by e2,max — (^e)^/4, where Ae is the specific la- 
tent heat. Figure ^ shows the behavior of 62 for various 
values of the system size. For L = 40 the curves deriving 
from histogram 2, from the combination of histograms 1 
and 3, and from 1 together with 3.2 to 3.5, respectively, 
all coincide, thereby confirming our approach. For the 
longer systems subphase histograms for the helical phase 
were used invariably. Included in the figure is the ex- 
trapolated function for infinite system size, correspond- 
ing to a latent heat of Ae = 1.37 Ehb and a transition 
at P = 2.78 Ehb- 

In closing, we will briefly discuss why other methods 
proposed for obtaining simulation results for first-order 
phase transitions fail in our system. Due to the high di- 
mensionality of the order parameter space and the com- 
plicated phase structure, compare Fig. |^, the straight- 
forward identification of a — possibly small number — 
of transition states in order parameter space, necessary 
for a successful application of the multicanonical method 
, is impossible. In addition, the states corresponding 
to the individual phases occupy only an extremely small 
part of the order parameter space. By raising the tem- 
perature, as in the entropic sampling Is) and simulated 
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tempering [|| methods, not only the transition states, 
but also an enormous number of "uninteresting" states 
would become accessible. To suppress these states at 
higher temperatures, a prohibitively large number of pa- 
rameters (of the order of 10^ for L — 100) would have 
to be optimized. The microcanonical ensemble approach 
promoted by Hiiller and coworkers (l^Jl^ or cluster dy- 
namics approaches [§,|l5| are also not applicable in our 
case. 

We believe that the approach presented here offers new 
possibilities for investigating model systems in parame- 
ter regimes where relaxation times are too long to allow 
equilibrated MC simulations, like at strong first-order 
transitions or at glass transitions and within the glass 
phase. A prerequisite for the successful application of 
the method is the existence of a detour in control pa- 
rameter space with a more gradual transition, i.e., much 
shorter relaxation times, so that the necessary free en- 
ergy differences can be determined. In order to apply 
the subphase method, Eqs. to (^, there has to be an 
order parameter set which allows to identify subphases 
unambiguously. 
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FIG. 1. Reconstructed histogram for L = 40, f3\EHB\ ~ 3.4 
(transition region), projected onto the (F,NhEi) plane. 

FIG. 2. Fluctuations in specific energy, 62 = ((e — (e))^^, 
vs P for L as indicated; also shown is the extrapolation for 
L —* 00. 



TABLE I. Histograms, L ^ 40. 



P Size 


Phase 


Correlation time 


[I-EhsT'] [lO^MC steps] 




[MC steps] 



Histograms, AS{h)/kB = -3.57 



1 


2.66 10 random coil 


30 


2 


3.43 600 mixed 


350000 


3 


3.96 50 multi-helical 


4600 




Auxiliary histograms, AS{h)/kB = —0.57 




la 


0.16 10 random coil 


50 


2a 


1.16 50 mixed 


200 


3a 


1.56 20 multi-helical 


300 



Subphase histograms, AS(/i)/fcs = —3.57 



3.2 3.96 


5 


Nh,i = 2 


100 


3.3 3.96 


5 


Nhel = 3 


100 


3.4 3.96 


5 


Nhel = 4 


160 


3.5 3.96 


5 


iVhei = 5 


170 


TABLE II. A/, L 


= 40, according to Eq. 


(0). 


Root mean square deviation 


(rmsd) according to Bennett Q. 


Histograms 




A/ 


rmsd 


1/2 




-1.627 


7-10-2 


2/3 




-22.884 


3 • 10-2 


E:l/3 




-24.511 


7.5 • 10-2 


1/la 




-9.501 


1.3 • 10-2 


la/2a 




-3.828 


8 • 10-^ 


2a/3a 




-14.741 


5.2 ■ 10-^ 


3a/3 




3.569 


9.3 ■ 10-^ 


E:l/3 




-24.501 


1.9 • 10-2 



